function [] = plot_SeasonalGrade(weekly_spl,time_spl,RWQname,rain,time_r)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   SEASONAL TIME-SERIES PLOT OF WEEKLY GRADE AND BOM RAINFALL DATA.     
%   SCATTER PLOT OF WEEKLY GRADE WITH COLOR CODED POINTS ACCORDING TO    
%   WATER QUALITY LEVEL.                                                     
%                                                                           
% INPUT     : - weekly_spl = weekly grade values for the season of interest          
%             - time_spl =  sampling date in matlab format
%             - RWQname = name of site (string) as in Bays_n_Beaches.mat
%             - rain = rainfall data over similar period as data             
%             - time_r = rainfall observation date   
%
% EXAMPLE   : load TestData.mat
%             load Bays_n_Beaches.mat
%             plot_SeasonalGrade(testdata.weekly5yr(:,5),...
%             testdata.timeweekly5yr(:,5),BEACHES{1},testdata.rain,...
%             testdata.rain_t)
%                                                                      
% Author: Benedicte Pasquer, IMOS/eMII (http://imos.org.au/)
% email: benedicte.pasquer@utas.edu.au
% May 2013   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin~=5
    error('Not enough input arguments')
end
if ~isnumeric(weekly_spl) ||~isnumeric(time_r)
    error('date must be in Matlab format')
end
if ~ischar(RWQname),error('RWQname must be a string'),end

% REMOVE MISSING DATA
weekly_spl(isnan(weekly_spl) | weekly_spl==99999) = [];
time_spl(isnan(time_spl) | time_spl==99999) = [];
% ADJUST PLOT XAXIS LIMITS FROM DEC 1 TIL MARCH 31 

date4period = datevec(time_spl);

dec1 = datenum(['1 dec ',num2str(date4period(1,1))]);
march31 = datenum(['31 march ', num2str(date4period(1,1)+1)]);
dectomarch = [dec1; time_spl; march31]; 

% GET RID OF RESAMPLING DATES FOR CLEAR DATE LABELING ON THE PLOT
for nd = 1:length(time_spl)
    tt = diff(time_spl);
    ind = find(tt<4);
    timeclean = time_spl ; 
    if ~isempty(ind)
        timeclean(ind+1)=[]; 
    end
end
    
% SORT INDEX OF DATA ACCORDING TO GRADE:
% GOOD : weekly_spl <140           
% RETEST: 140< weekly_spl <280
% FAIL : weekly_spl >280
  
goodrwq = (weekly_spl <140);
retestrwq = (weekly_spl>=140 & weekly_spl <280);
failrwq = (weekly_spl >280);
    
% ARRAY OF DATA TO DISPLAY RESAMPLING AND POTENTIAL BEACH CLOSURE LEVEL
GoodQ = ones(1,length(dectomarch)) * 140;    
BadQ = ones(1,length(dectomarch)) * 280;     
    
% FIGURE

figure('position',[1 1 1500 900],'Color',[.9 .98 .98])
hold on

l1 = line(dectomarch,GoodQ,'Linestyle','--','Linewidth',1);
l2 = line(dectomarch,BadQ,'Linestyle','-.','Linewidth',1);

% 'GOOD' RESULTS DISPLAY IN GREEN 
if  ~isempty(find(goodrwq==1)) 
    l3 = line(time_spl(goodrwq),weekly_spl(goodrwq),...
        'Linestyle','none','Marker','o',...
        'MarkerSize',10,'MarkerFaceColor','g','MarkerEdgeColor','k');
else 
    l3=[];
end
% 'RETEST' RESULTS DISPLAY IN YELLOW
if ~isempty(find(retestrwq==1))
    l4 = line(time_spl(retestrwq),weekly_spl(retestrwq),...
        'Linestyle','none','Marker','o',...
        'MarkerSize',10,'MarkerFaceColor','y','MarkerEdgeColor','k');
else 

    l4=[];
end
% 'FAIL' RESULTS DISPLAY IN RED 
if ~isempty(find(failrwq==1))
    l5 = line(time_spl(failrwq),weekly_spl(failrwq),...
        'Linestyle','none','Marker','o',...
        'MarkerSize',10,'MarkerFaceColor','r','MarkerEdgeColor','k');
else
    l5 = [];
end
     
%TITLE AND LABELS
titre{1} = [num2str(date4period(1,1)) '-' num2str(date4period(1,1)+1) ' water quality trend'];
titre{2} = ['at ', RWQname]; 
titre{3} = ' ';

title(titre,...
    'Fontname','arial','Fontsize',14,'Fontweight','bold', ...
    'color',[0 0 .7]);

set(gca, 'Xlim', [dec1-2 march31], 'Xtick',timeclean,...
    'Color', [.9 .98 .98],'YColor', [0 0 .7], 'Xcolor', [0 0 .7], ...       
    'Fontname','arial','Fontsize',14,'Fontweight','bold',...
    'Position' , [0.27 0.2 0.675 0.70],'TickDir','out')
datetick('x','dd/mm/yyyy','keepticks')
hh = rotateXLabels(gca(),45);
set (hh,'Color',[0 0 .7])
xlabel('Sampling date')
ylabel('Bacterial level [Enterococci/100ml]')
set(gca, 'XLim',[dec1-2 march31])
hold off

set(gcf, 'InvertHardCopy', 'off');

% 
%ADD A SECOND AXE FOR BOM RAINFALL DATA
%  
    
ax1 = gca;
ax2 = axes('Position',get(ax1,'Position'),...
    'XAxisLocation','top',...
    'YAxisLocation','right',...
    'Color','none',...
    'XColor','k','YColor','k',...
    'XTick',[ ],...
    'XLim',[dec1-2 march31],... 
    'YLim',[0 max(rain)*1.25], ... %may need to adjust YLim based on rain data value
    'Fontname','arial','Fontsize',14,'Fontweight','bold');

l6 = line(time_r,rain,...
    'Color','k','Marker','d','Parent',ax2);
ylabel(ax2, 'Rainfall [mm]')
 
%DISPLAY LEGEND OF PLOTTED DATA ONLY
% list of line handle   
line_h = [l1 l2 l3 l4 l5 l6];

Dyn = {char('Trigger level','for resampling'),...
    char('Trigger level','for potential beach closure'),...
    'Good','Retest','Fail',...
    char('Rainfall (Hobart, Ellerslie Rd)','source: Bureau of Meteorology')};

in=[];
for j = 1:6
    if ~eval(['isempty(l', num2str(j),')'])
        in = [in ;j];
    end
end

hleg = legend(line_h,Dyn(in),'Location', [.07 .6 .07 .2]);
set(hleg,'Fontsize',12)

% OUTPUT FILE
figname = strcat(RWQname,'ShortTermGrade');
figname(isspace(figname)) = []; %remove blank from filename
saveas(gcf,figname,'png')




